// <copyright file="SpecialFunctionsTests.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2010 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>

namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
{
    using System;
    using NUnit.Framework;

    /// <summary>
    /// Special functions tests.
    /// </summary>
    [TestFixture, Category("Functions")]
    public class SpecialFunctionsTests
    {
        /// <summary>
        /// Di-Gamma function
        /// </summary>
        /// <param name="x">Input X value.</param>
        /// <param name="f">Function result.</param>
        [TestCase(Double.NaN, Double.NaN)]
        [TestCase(-1.5, 0.70315664064524318722569033366791109947350706200623256)]
        [TestCase(-0.5, 0.036489973978576520559023667001244432806840395339565891)]
        [TestCase(0.1, -10.423754940411076232100295314502760886768558023951363)]
        [TestCase(1.0, -0.57721566490153286060651209008240243104215933593992359)]
        [TestCase(1.5, 0.036489973978576520559023667001244432806840395339565888)]
        [TestCase(Constants.Pi / 2, 0.10067337642740238636795561404029690452798358068944001)]
        [TestCase(2.0, 0.42278433509846713939348790991759756895784066406007641)]
        [TestCase(2.5, 0.70315664064524318722569033366791109947350706200623255)]
        [TestCase(3.0, 0.92278433509846713939348790991759756895784066406007641)]
        [TestCase(Constants.Pi, 0.97721330794200673329206948640618234364083460999432603)]
        [TestCase(3.5, 1.1031566406452431872256903336679110994735070620062326)]
        [TestCase(4.0, 1.2561176684318004727268212432509309022911739973934097)]
        [TestCase(4.5, 1.3888709263595289015114046193821968137592213477205183)]
        [TestCase(5.0, 1.5061176684318004727268212432509309022911739973934097)]
        [TestCase(5.5, 1.6110931485817511237336268416044190359814435699427405)]
        [TestCase(10.1, 2.2622143570941481235561593642219403924532310597356171)]
        public void DiGamma(double x, double f)
        {
            AssertHelpers.AlmostEqualRelative(f, SpecialFunctions.DiGamma(x), 12);
        }

        /// <summary>
        /// Di-Gamma inverse.
        /// </summary>
        /// <param name="x">Input X value.</param>
        /// <param name="f">Function result.</param>
        [TestCase(Double.NaN, Double.NaN)]
        [TestCase(0.0, Double.NegativeInfinity)]
        [TestCase(0.1, -10.423754940411076232100295314502760886768558023951363)]
        [TestCase(1.0, -0.57721566490153286060651209008240243104215933593992359)]
        [TestCase(1.5, 0.036489973978576520559023667001244432806840395339565888)]
        [TestCase(Constants.Pi / 2, 0.10067337642740238636795561404029690452798358068944001)]
        [TestCase(2.0, 0.42278433509846713939348790991759756895784066406007641)]
        [TestCase(2.5, 0.70315664064524318722569033366791109947350706200623255)]
        [TestCase(3.0, 0.92278433509846713939348790991759756895784066406007641)]
        [TestCase(Constants.Pi, 0.97721330794200673329206948640618234364083460999432603)]
        [TestCase(3.5, 1.1031566406452431872256903336679110994735070620062326)]
        [TestCase(4.0, 1.2561176684318004727268212432509309022911739973934097)]
        [TestCase(4.5, 1.3888709263595289015114046193821968137592213477205183)]
        [TestCase(5.0, 1.5061176684318004727268212432509309022911739973934097)]
        [TestCase(5.5, 1.6110931485817511237336268416044190359814435699427405)]
        [TestCase(10.1, 2.2622143570941481235561593642219403924532310597356171)]
        public void DiGammaInv(double x, double f)
        {
            AssertHelpers.AlmostEqualRelative(x, SpecialFunctions.DiGammaInv(f), 13);
        }

        /// <summary>
        /// Compute the t'th harmonic number using a loop.
        /// </summary>
        /// <param name="t">Input value.</param>
        /// <returns>Harmonic number.</returns>
        private static double ExactHarmonic(int t)
        {
            var r = 0.0;
            for (var i = 1; i <= t; i++)
            {
                r += 1.0 / i;
            }

            return r;
        }

        /// <summary>
        /// Harmonic function.
        /// </summary>
        /// <param name="i">Input value.</param>
        [TestCase(1)]
        [TestCase(2)]
        [TestCase(4)]
        [TestCase(8)]
        [TestCase(16)]
        [TestCase(100)]
        [TestCase(1000)]
        [TestCase(10000)]
        [TestCase(100000)]
        [TestCase(1000000)]
        public void Harmonic(int i)
        {
            AssertHelpers.AlmostEqualRelative(ExactHarmonic(i), SpecialFunctions.Harmonic(i), 12);
        }

        /// <summary>
        /// Log beta function.
        /// </summary>
        [Test]
        public void BetaLn()
        {
            AssertHelpers.AlmostEqualRelative(Math.Log(0.5), SpecialFunctions.BetaLn(1.0, 2.0), 13);
            AssertHelpers.AlmostEqualRelative(Math.Log(1.0), SpecialFunctions.BetaLn(1.0, 1.0), 13);
        }

        /// <summary>
        /// Beta function.
        /// </summary>
        [Test]
        public void Beta()
        {
            AssertHelpers.AlmostEqualRelative(0.5, SpecialFunctions.Beta(1.0, 2.0), 13);
            AssertHelpers.AlmostEqualRelative(1.0, SpecialFunctions.Beta(1.0, 1.0), 13);
        }

        /// <summary>
        /// Beta incomplete.
        /// </summary>
        /// <param name="a">A parameter.</param>
        /// <param name="b">B parameter.</param>
        /// <param name="x">Input X value.</param>
        /// <param name="f">Function result.</param>
        [TestCase(0.100000, 0.100000, 0.100000, 8.0117356206774655704238957309013421730449536344797)]
        [TestCase(0.100000, 0.100000, 0.500000, 9.8573197445250802696495512442154091185464210677262)]
        [TestCase(0.100000, 0.100000, 0.800000, 11.045931323774722512127911008108711428206559153167)]
        [TestCase(0.100000, 1.500000, 0.100000, 7.906686887040059574762793470136627722332467302241)]
        [TestCase(0.100000, 1.500000, 0.500000, 9.1012542128394916083902345366778109992938115490192)]
        [TestCase(0.100000, 1.500000, 0.800000, 9.368806890135865087467208304972858639790010273545)]
        [TestCase(0.100000, 2.500000, 0.100000, 7.8363997919172974609110616787835466667385876309684)]
        [TestCase(0.100000, 2.500000, 0.500000, 8.7385989355952880797816267602947881842609716646237)]
        [TestCase(0.100000, 2.500000, 0.800000, 8.8379245631995373736838511405978394101548418300997)]
        [TestCase(0.100000, 5.500000, 0.100000, 7.6464829466025722609741358285691936422732595850987)]
        [TestCase(0.100000, 5.500000, 0.500000, 8.0832162883698521112308593494639236817925940654185)]
        [TestCase(0.100000, 5.500000, 0.800000, 8.089843275520988350320569916795581670946930566783)]
        [TestCase(1.500000, 0.100000, 0.100000, 0.022303834332028031481145155068151014763064226606192)]
        [TestCase(1.500000, 0.100000, 0.500000, 0.33465159984030260689318592198561677643426763544959)]
        [TestCase(1.500000, 0.100000, 0.800000, 1.002080089012839104050394134134167486035896540574)]
        [TestCase(1.500000, 1.500000, 0.100000, 0.02043763859916055001568569052740016093388906615616)]
        [TestCase(1.500000, 1.500000, 0.500000, 0.19634954084936207740391521145496893026232308746094)]
        [TestCase(1.500000, 1.500000, 0.800000, 0.33678717944852264351783475904713816723851995610486)]
        [TestCase(1.500000, 2.500000, 0.100000, 0.019218819299580275673976660038794008316192763455412)]
        [TestCase(1.500000, 2.500000, 0.500000, 0.13984143709134770536862427239415113179782821039714)]
        [TestCase(1.500000, 2.500000, 0.800000, 0.18972692305759464976318019465615085035583828745353)]
        [TestCase(1.500000, 5.500000, 0.100000, 0.016056550082674778556355475820652277420877772570937)]
        [TestCase(1.500000, 5.500000, 0.500000, 0.061380263212265132490746506045997506787829048203228)]
        [TestCase(1.500000, 5.500000, 0.800000, 0.064403479961606576649564368278872635247235510624673)]
        [TestCase(2.500000, 0.100000, 0.100000, 0.001352753157951915161848451666929681258878645850634)]
        [TestCase(2.500000, 0.100000, 0.500000, 0.10756276379201896635529117527407904299772826375257)]
        [TestCase(2.500000, 0.100000, 0.800000, 0.5587192957063609402827988523062144494394362102048)]
        [TestCase(2.500000, 1.500000, 0.100000, 0.0012188192995802743417090304886061526176963027007477)]
        [TestCase(2.500000, 1.500000, 0.500000, 0.056508103758014372035290939060817798464494877063805)]
        [TestCase(2.500000, 1.500000, 0.800000, 0.14706025639092799375465456439098731688268166865133)]
        [TestCase(2.500000, 2.500000, 0.100000, 0.0011320572373426029655709496224583878329664195067652)]
        [TestCase(2.500000, 2.500000, 0.500000, 0.036815538909255389513234102147806674424185578898927)]
        [TestCase(2.500000, 2.500000, 0.800000, 0.067947596146597995171095886486269312250373204235372)]
        [TestCase(2.500000, 5.500000, 0.100000, 0.00091001787485888099434748885472859856478675247779447)]
        [TestCase(2.500000, 5.500000, 0.500000, 0.012036842116913956962302822724142322883106224614977)]
        [TestCase(2.500000, 5.500000, 0.800000, 0.0137861171346299807272679372975664758815098236357)]
        [TestCase(5.500000, 0.100000, 0.100000, 0.00000062268687636453588281206442354375188797585844216401)]
        [TestCase(5.500000, 0.100000, 0.500000, 0.0066577573700032687517874433540471831930305154614716)]
        [TestCase(5.500000, 0.100000, 0.800000, 0.15893826959760073090597940189540380779121963557728)]
        [TestCase(5.500000, 1.500000, 0.100000, 0.00000055008267477746389601958949824166456746390970113234)]
        [TestCase(5.500000, 1.500000, 0.500000, 0.0030469298789317991574131727126641734544957148698945)]
        [TestCase(5.500000, 1.500000, 0.800000, 0.029928813294939917230433606365228383142081556070117)]
        [TestCase(5.500000, 2.500000, 0.100000, 0.00000050358914459517094905570885392504304689450166185919)]
        [TestCase(5.500000, 2.500000, 0.500000, 0.0017689849740568141051599655812851800259633674721202)]
        [TestCase(5.500000, 2.500000, 0.800000, 0.010158231420344267874007806875642686140428489302542)]
        [TestCase(5.500000, 5.500000, 0.100000, 0.00000038687628134504851234251462244340434107233073666673)]
        [TestCase(5.500000, 5.500000, 0.500000, 0.00037750308451873202137593561772653328266987165863158)]
        [TestCase(5.500000, 5.500000, 0.800000, 0.00074361660080007708142054771607396897718180545812717)]
        public void BetaIncomplete(double a, double b, double x, double f)
        {
            AssertHelpers.AlmostEqualRelative(f, SpecialFunctions.BetaIncomplete(a, b, x), 11);
        }

        /// <summary>
        /// Beta regularized.
        /// </summary>
        /// <param name="a">A parameter.</param>
        /// <param name="b">B parameter.</param>
        /// <param name="x">Input X value.</param>
        /// <param name="f">Function result.</param>
        [TestCase(0.100000, 0.100000, 0.100000, 0.40638509393627598963947434031370208398700911383034)]
        [TestCase(0.100000, 0.100000, 0.500000, 0.5)]
        [TestCase(0.100000, 0.100000, 0.800000, 0.56029080977665439586092261707431380702999082404228)]
        [TestCase(0.100000, 1.500000, 0.100000, 0.83793618164513694339361766577607270440777026543289)]
        [TestCase(0.100000, 1.500000, 0.500000, 0.96453423693668031005965611526662377555163072304903)]
        [TestCase(0.100000, 1.500000, 0.800000, 0.99288897919542998751072045939791723468184549496473)]
        [TestCase(0.100000, 2.500000, 0.100000, 0.88585310309894667823419197676046986281131524104959)]
        [TestCase(0.100000, 2.500000, 0.500000, 0.98784074184406228317625751436240720987391458581291)]
        [TestCase(0.100000, 2.500000, 0.800000, 0.99906884630106408932550166100889714490764416359312)]
        [TestCase(0.100000, 5.500000, 0.100000, 0.94519184147610137383076366494492317664691771837548)]
        [TestCase(0.100000, 5.500000, 0.500000, 0.99917702583101287494949180374093750006040749104055)]
        [TestCase(0.100000, 5.500000, 0.800000, 0.99999619645266501604882303584892931538704382134298)]
        [TestCase(1.500000, 0.100000, 0.100000, 0.0023637194748231330327322712604688049136736986021871)]
        [TestCase(1.500000, 0.100000, 0.500000, 0.03546576306331968994034388473337622444836927695097)]
        [TestCase(1.500000, 0.100000, 0.800000, 0.10619861080705813882414161079467267590245761692272)]
        [TestCase(1.500000, 1.500000, 0.100000, 0.052044019330913933551809009997594089261007490064035)]
        [TestCase(1.500000, 1.500000, 0.500000, 0.5)]
        [TestCase(1.500000, 1.500000, 0.800000, 0.85762151006735301922188173009147886735366387396138)]
        [TestCase(1.500000, 2.500000, 0.100000, 0.097880642941379793645839194076629344959184516635409)]
        [TestCase(1.500000, 2.500000, 0.500000, 0.71220659078919378102517835116335248271261286098728)]
        [TestCase(1.500000, 2.500000, 0.800000, 0.96627128455142020796603976406510565066301101575357)]
        [TestCase(1.500000, 5.500000, 0.100000, 0.2492200779249636429835386068413815093625886315581)]
        [TestCase(1.500000, 5.500000, 0.500000, 0.95270739368361339952038048248181862978690743677285)]
        [TestCase(1.500000, 5.500000, 0.800000, 0.99963193911681378117173263163355516583213310848701)]
        [TestCase(2.500000, 0.100000, 0.100000, 0.00015291978644767572330394885397661411972523657814463)]
        [TestCase(2.500000, 0.100000, 0.500000, 0.012159258155937716823742485637592790126085414187091)]
        [TestCase(2.500000, 0.100000, 0.800000, 0.063159516487818477363383984899205310365144861115492)]
        [TestCase(2.500000, 1.500000, 0.100000, 0.006207395720448073457778825918558833562830463492662)]
        [TestCase(2.500000, 1.500000, 0.500000, 0.28779340921080621897482164883664751728738713901272)]
        [TestCase(2.500000, 1.500000, 0.800000, 0.74897173558328583047772369611785208404431673216919)]
        [TestCase(2.500000, 2.500000, 0.100000, 0.015374720442541245985473611768528587536206924004576)]
        [TestCase(2.500000, 2.500000, 0.500000, 0.5)]
        [TestCase(2.500000, 2.500000, 0.800000, 0.92281137475779334211841505067903343661808086762039)]
        [TestCase(2.500000, 5.500000, 0.100000, 0.065915491253258347344782598110283193803051483319306)]
        [TestCase(2.500000, 5.500000, 0.500000, 0.87186678766868243532031253918149387446781682306342)]
        [TestCase(2.500000, 5.500000, 0.800000, 0.99857234512565487924356478796968439529913692090276)]
        [TestCase(5.500000, 0.100000, 0.100000, 0.000000076971146008440526509370378069895435149621047562911)]
        [TestCase(5.500000, 0.100000, 0.500000, 0.00082297416898712505050819625906249993959250895944939)]
        [TestCase(5.500000, 0.100000, 0.800000, 0.01964656911825443767408345737200950450690987324177)]
        [TestCase(5.500000, 1.500000, 0.100000, 0.0000085380512231662773545327104780279605728164994539023)]
        [TestCase(5.500000, 1.500000, 0.500000, 0.047292606316386600479619517518181370213092563227146)]
        [TestCase(5.500000, 1.500000, 0.800000, 0.46453697358156781101222907948783656601193056166297)]
        [TestCase(5.500000, 2.500000, 0.100000, 0.000036476564661926426853537763285519105619982887967225)]
        [TestCase(5.500000, 2.500000, 0.500000, 0.12813321233131756467968746081850612553218317693658)]
        [TestCase(5.500000, 2.500000, 0.800000, 0.73579303531824700577792236664278017906658620718916)]
        [TestCase(5.500000, 5.500000, 0.100000, 0.00051241472879389331855418481833750204851457778933154)]
        [TestCase(5.500000, 5.500000, 0.500000, 0.5)]
        [TestCase(5.500000, 5.500000, 0.800000, 0.98491460241721309518021565227501009891015897915977)]
        public void BetaRegularized(double a, double b, double x, double f)
        {
            AssertHelpers.AlmostEqualRelative(f, SpecialFunctions.BetaRegularized(a, b, x), 11);
        }

        [Test]
        public void BetaRegularizedLargeArguments()
        {
            Assert.That(
                SpecialFunctions.BetaRegularized(3846.2382301675848, 738420369.64263761, 0.0000052331266889206809),
                Is.EqualTo(0.61624368331298802128).Within(1e-5));
        }

        /// <summary>
        /// Logit function.
        /// </summary>
        /// <param name="p">Input value.</param>
        /// <param name="x">Expected X value.</param>
        [TestCase(0.000000, Double.NegativeInfinity)]
        [TestCase(0.000010, -11.512915464920228103874353849992239636376994324587)]
        [TestCase(0.001000, -6.9067547786485535272274487616830597875179908939086)]
        [TestCase(0.100000, -2.1972245773362193134015514347727700402304323440139)]
        [TestCase(0.500000, 0.0)]
        [TestCase(0.900000, 2.1972245773362195801634726294284168954491240598975)]
        [TestCase(0.999000, 6.9067547786485526081487245019905638981131702804661)]
        [TestCase(0.999990, 11.512915464924779098232747799811946290419057060965)]
        [TestCase(1.000000, Double.PositiveInfinity)]
        public void Logit(double p, double x)
        {
            AssertHelpers.AlmostEqualRelative(x, SpecialFunctions.Logit(p), 15);
        }

        /// <summary>
        /// Logistic function.
        /// </summary>
        /// <param name="p">Expected value.</param>
        /// <param name="x">Input X value.</param>
        [TestCase(Double.NegativeInfinity, 0.000000)]
        [TestCase(-11.512915464920228103874353849992239636376994324587, 0.000010)]
        [TestCase(-6.9067547786485535272274487616830597875179908939086, 0.001000)]
        [TestCase(-2.1972245773362193134015514347727700402304323440139, 0.100000)]
        [TestCase(0.0, 0.500000)]
        [TestCase(2.1972245773362195801634726294284168954491240598975, 0.900000)]
        [TestCase(6.9067547786485526081487245019905638981131702804661, 0.999000)]
        [TestCase(11.512915464924779098232747799811946290419057060965, 0.999990)]
        [TestCase(Double.PositiveInfinity, 1.000000)]
        public void Logistic(double p, double x)
        {
            AssertHelpers.AlmostEqualRelative(x, SpecialFunctions.Logistic(p), 14);
        }
    }
}
